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Abstract 

Residual distributions (RD) schemes are a class of of high-resolution finite volume 
methods for unstructured grids. A key feature of these schemes is that they make 
use of genuinely multidimensional (approximate) Riemann solvers as opposed to 
the piecemeal ID Riemann solvers usually employed by finite volume methods. In 
ID, LeVeque and Pelanti [J. Com,p. Phys. 172, 572 (2001)] showed that many of 
the standard approximate Riemann solver methods (e.g., the Roe solver, HLL, Lax- 
Priedrichs) can be obtained from applying an exact Riemann solver to relaxation 
systems of the type introduced by Jin and Xin [Comm. Pure Appl. Math. 48, 235 
(1995)]. In this work we extend LcVcquc and Polanti's results and obtain a multi- 
dimensional relaxation system from which multidimensional approximate Riemann 
solvers can be obtained. In particular, we show that with one choice of parameters 
the relaxation system yields the standard N-scheme. With another choice, the re- 
laxation system yields a new Riemann solver, which can be viewed as a genuinely 
multidimensional extension of the local Lax-Friedrichs scheme. This new Riemann 
solver does not require the use Roe-Struijs-Deconinck averages, nor does it require 
the inversion of an m x m matrix in each computational grid cell, where m is the 
number of conserved variables. Once this new scheme is established, we apply it on 
a few standard cases for the 2D compressible Euler equations of gas dynamics. We 
show that through the use of linear-preserving limiters, the new approach produces 
numerical solutions that are comparable in accuracy to the N-scheme, despite being 
computationally less expensive. 
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1 Introduction 



In the last few decades intense researcli into sliock-capturing schemes has 
resulted in several numerical methods for solving partial differential equa- 
tions (PDEs) that admit discontinuous weak solutions in the form of shock- 
waves. Examples of such schemes include WENO (weighted essentially non- 
oscillatory) [TB], central [12], MUSCL (monotone upstream-centered schemes 
for conservation laws) [27] , and wave propagation schemes [I9] . One difficulty 
with these methods is that in general they do not trivially extend to problems 
in complex geometries. In order to handle application problems where com- 
plex geometry is of great importance, three broad classes of strategies have 
been considered: (1) Cartesian cut-cell methods jT3], (2) overlapping meshes 
in] , and (3) unstructured grid methods. We will focus in this work on the last 
approach, thus eliminating the problem of cut-cells (1st approach) and inter- 
polation between different grid patches (2nd approach), but requiring some 
efficient grid generation tool. 

On unstructured grids, the two main classes of methods that have been devel- 
oped are discontinuous Galerkin (DG) [T0f22] and residual distribution (RD) 
[Tp3] schemes. The discontinuous Galerkin approach is based on defining a 
piecewise polynomial approximation that is continuous inside element interi- 
ors, but discontinuous across element boundaries. Local ID Riemann problems 
are solved across element boundaries to construct the necessary numerical 
fluxes. Residual distribution schemes can be viewed as a finite volume method 
where the finite volumes are defined by a grid that is dual to the original 
triangulation. 

Although DG schemes have their own particular advantages, the focus of this 
work will be on RD schemes and, in particular, the aspect of RD schemes 
that separates them from all other methods: RD schemes are based on solving 
genuinely multi-dimensional Riemann problems. This aspect allows one to 
obtain methods that are positivity preserving for scalar conservation laws and 
essentially non-oscillatory for systems. This same feature, however, presents 
a challenge: how can these multi-dimensional Riemann problems be solved 
efficiently? The standard answer to this question is the so-called systems N- 
scheme |26] (see also [I]|3|), which is a generalization of Roe's approximate 
Riemann solver for ID systems [23]. One goal of this work is to develop an 
alternative to this approach. 

LeVeque and Pelanti [21] showed how several of the standard approximate 
Riemann solvers can be interpreted as exact Riemann solvers for a perturbed 
system of hyperbolic equations known as relaxation systems. Their work was 
motivated by Jin and Xin's earlier paper ^17i| on a class of numerical methods 
known as relaxation schemes. What LeVeque and Pelanti essentially showed 
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is that Jin and Xin's "new" class of methods could actually be thought of as 
a reinterpretation of various pre-existing approximate Riemann solvers; these 
results are reviewed in Section [2j After reviewing RD schemes in Section |3} 
we focus in this work on the continuation of LeVeque and Pelanti's reason- 
ing and show how the N-scheme can be also be derived from a relaxation 
system. Furthermore, using this interpretation we derive a novel genuinely 
multi-dimensional Riemann solver that can be viewed as a multidimensional 
extension of the ID local Lax-Friedrichs scheme [24j. Both of these results are 
presented in Section |4j Finally, we compare the numerical accuracy of the N- 
scheme and the newly derived scheme on several examples in Section |5| What 
we find is that when the appropriate limiters are applied, the novel scheme has 
comparable accuracy to the N-scheme, although it tends to be slightly more 
diffusive - this result is of course consistent with well-known ID results com- 
paring local Lax-Friedrichs versus Roe-type approximate Riemann solvers. On 
the other hand, this loss of accuracy is compensated by the fact that the new 
scheme is less computational expensive. This gain in computational efficiency 
will become significant for problems involving complicated equations such as 
the relativistic Euler or MHD equations. 



2 Review of ID relaxation systems 

We briefiy review in this section the results of LeVeque and Pelanti [21] for 
the case of ID conservation laws. For simplicity we consider for the moment 
a scalar conservation laws of the form 

q,t + f{q),. = 0, (1) 

where x G M is the spatial coordinate, t G M"*" is the time coordinate, g G M 
is the conserved variable, and f{q) : M M is the fiux function. We assume 
that this conservation law is hyperbolic, meaning that /'(<?) G M for all q in 
the solution domain. 

2.1 Finite volume methods in ID 

Using the idea of relaxation, we will construct in this section numerical meth- 
ods for approximating ([T]). All of these methods are in the general class of 
finite volume methods [20], which we briefiy recall in this subsection. 

Let Tax be the numerical grid with grid cells centered at x = Xi and spanning 
the interval [xi — Ax/2,Xi + Ax/2], where 

Xi = a + {i -1/2) Ax. (2) 
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Here i is an integer ranging from 1 to A^, a and b are the left and right end 
points of the domain, respectively, and Ax = (6 — a)/A^ is the grid spacing. 
In each grid cell Xi and at each time level t = t" we seek an approximation to 
the cell average of the exact solution q{x,t): 

1 i'Xi+Ax/2 
Ax Jxi-Ax/2 

Integrating ([T]) over the grid cell centered at Xj and from t = to t = 
t"+i results in a numerical update formula for that can be written in the 
following fluctuation splitting form: 

=Q7-^ \A-AQI\,,, + A^AQUJ , (4) 



Ax 



where A~ AQ^_^_^^2 •A^AQ^_^^2 l^ft- and right-going fluctuations, which 
measure the amount of flux that enters into grid cell Xj through the grid 
interfaces at x = Xj + Ax/2 and x = Xj — Ax/2, respectively. In order for this 
update to be numerically conservative these fluctuations must satisfy 

A-AQ-^,/2 + ^+AQIVi/2 = f{Q7+i) - fm- (5) 

Note that in update Q we collect the left-going fluctuation from the grid 
interface at Xi + Ax/2 and the right-going fluctuation from the grid interface 
at Xi — Ax/2, while in expression ([s]) we are adding the left- and right-going 
fluctuations at the same grid interface. 

For first-order accurate methods, the fluctuations in update Q are obtained 
by first assuming that the approximate solution has a constant value, Q", in 
each grid cell, and then solving at each grid interface, Xj_i/2 = Xi — Ax/2, the 
initial value problem for ([T]) with the piecewise constant initial data: 

[Q" if X < Xi_i/2. 

This initial value problem is referred to as the Riemann problem. One of the 
pieces of information that can be obtained from solving the Riemann problem 
is how much of the initial flux difference, /(Q") — is carried to the 

left and how much to the right. It is precisely this information that is stored 
in the fluctuations, A~AQ and A'^AQ. 



2.2 Relaxation method framework in ID 



The relaxation schemes introduced by Jin and Xin are based on the idea of 
approximating the quasilinear equation ([T]) by a linear system with a cleverly 
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chosen source term. The role of this source term is to force the hnear system 
to relax in the hmit as t — >■ oo towards the original equation. By "hiding" the 
nonlinearity in the source term, relatively complicated quasilinear Riemann 
problems can be replaced by simpler linear Riemann problems. 



There are many kinds of relaxation systems that one could develop in order 
to create an approximate solution to ([T]) (see pp. 26-48 of Bouchut [7] for a 
discussion of several different approaches). In this work we follow the approach 
of [21] and consider the following relaxation system: 
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source term 



where c, c? G M are parameters that will be adjusted in the next few subsections 
in order to arrive at various approximate Riemann solvers. Without loss of 
generality we will assume that c < d. The key observation is that by taking 
e — >• 0, the right-hand side forces fi /(<?)• Since the first equation in the 
above system is g^j + = 0, /i /(g) will cause the relaxed system solution 
to approach the original conservation law solution. 

In order to make this statement more precise, we will carry out a so-called 
Chapman- Enskog expansion, which in this case is simply a Taylor series ex- 
pansion in € applied to system ([t]). Omitting the algebra, this expansion to 
0{e'^) yields the following equation for q{x,t): 

q,t + f{q),x =£ 
'■ V ' 

original cons, law s,. 

This approximation is stable for e > if the values c and d are chosen to 
produce positive (or at least non-negative) diffusion; this occurs if 

c<|<<;. (9) 

The above statement is often referred as the suh- characteristic condition (see 
for example Pll7|21] ) , since it requires that the eigenvalues of the coefficient 
matrix, which are just c and d, enclose the characteristic speed of the original 
conservation law, df jdq. 

From relaxation system ([T]), LeVeque and Pelanti [21] showed that various 
classical approximate Riemann solvers could be derived. Following the phi- 
losophy of operator splitting (see pp. 380 — 390 of LeVeque [20] for a review). 
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system ([T]) is first rewritten as two sub-problems: 
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(10) 

(11) 



Using this interpretation, LeVeque and Pelanti's |21j procedure for obtaining 
different approximate Riemann solvers can be summarized as follows: 

(1) Choose values for the parameters c and d. 

(2) Exactly solve the Riemann problem for the homogeneous linear system 



(10). 



(3) Approximate the effect of equation ( 11 ) on the solution calculated in Step 
(2) by directly setting /i = /(g). In other words, instantaneously relax 
the solution from Step (2) to the e — limit. 

We will simply refer to this as the relaxation procedure. In the next four sub- 
sections, we will apply this strategy for various values of c and d. Each time 
we carry out step (2) of the above procedure we will exactly solve the initial 
value problem (i.e., — oo < a; < oo) for system (10) using the generic Riemann 
data: 



q{x,0) 



if X < 0, 
if X > 0, 



and /i(x, 0) 




if a; < 0, 
if a; > 0, 



(12) 



where Qi and Qr are constants. Note that we are allowed to take = f{q) in 
the initial conditions at any arbitrary time step, since in the previous time-step 
we set = /(g) in Step (3) of the relaxation procedure. 



2.3 Local Lax-Friedrichs (LLF) for scalar equations 



The local Lax-Friedrichs or Rusanov method ^24j is obtained by applying the 
relaxation procedure with the choice 

c = —s and d = s, (13) 

where s > |/'(g)| in order to satisfy the sub-characteristic condition. With 
this choice the Riemann solution is obtained by splitting the jump between 
the left and right states, {Qe, f{Qe)) and (Qr, fiQr)), along the eigenvectors 
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of the coefficient matrix: 



Qr 


-Qi 
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(14) 



where the corresponding eigenvalues are = — s and = s. From this 
expression we obtain the following fluctuations: 



A'AQ = X^a^ 



2 ifiQr) - fiQi)) - I (Qr - Qi) 



1 



ifiQr) - fiQi)) + ^iQr- Qi) . 



(15) 
(16) 



2.4 Harten, Lax, and van Leer (HLL) for scalar equations 



The HLL method of [H] is obtained by applying the above procedure with 
the choice 



c = se and d = Sr 



(17) 



where < f'iq) < Sr in order to satisfy the sub-characteristic condition. With 
this choice the Riemann solution is obtained by splitting the jump between 
the left and right states, iQi, fiQe)) and iQr, fiQr)), along the eigenvectors 
of the coefficient matrix: 



Qr 


-Qi 
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+ 


1 


fiQr) 


- fiQi) 




Si 







where the corresponding eigenvalues are A^ = se and A^ = Sr- From this 
expression we obtain the following fluctuations: 



ifiQr 

- Si I 



- fiQi)) 
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Sn s 



i -'r 



iQr — Qi 



(19) 



V Sr - Si 

In the above expressions we have made use of the following notation: 

s+=max(0,s) and s~=min(0, s). (20) 
We will make use of this notation throughout the remainder of this paper. 
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2. 5 Roe 's approximate Riemann solver for scalar equations 



Roe's approximate Riemann solver 
procedure with the choice 

c = d = s = 



is obtained by applying the above 

fm - fiQi) 



Qr ~ Qe 



(21) 



With this choice the coefficient matrix becomes deficient, since only one lin- 
early independent eigenvector exists. Therefore, the jump between the left and 
right states, {Qi,f{Qe)) and {Qr, f{Qr)), can be written as 



(22) 



where the corresponding eigenvalue is A = s (algebraic multiplicity 2, geomet- 
ric multiplicity 1). Although this seems like an over-determined system for 
a, there exists a unique solution: a = Qr — Qe- This results in the following 
fiuctuations: 



Qr 


~Qi 


= a 


1 


fm 


- fiQi) 




s 



A^AQ = {Qr - Qi 



(23) 



2. 6 LLF and HLL for systems 



Finally, we briefiy explain how these three interpretations can be applied to a 
system of conservation laws of the form ([T|, now with q G and f{q) : — 
M*". We again assume hyperbolicity, which implies that the m x m matrix, 
df/dq, has m real eigenvalues and m linearly independent eigenvectors for all 
q in the solution domain. 

The systems LLF and HLL methods are obtained by considering the following 
relaxation system: 
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cdl {c + d)I 
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(24) 



where I is the m x m identity matrix, 01 is the m x m matrix with zeros in 
every entry, /i G M™, and c, (i G M. 

The systems LLF method is obtained by taking 

s = d = —c, where s> max |A^|, (25) 

p=l,...,m 
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and A'' is the p^^ eigenvalue of df/dq. With this choice we again arrive at 
formula (15), which is now applied to each component of the solution vector. 

Similarly, the systems HLL method is obtained by taking 

S£ = c, Sr = d, where S£ < min (A*^) and Sr > max (A^) . (26) 

p=l,...,m p=l,...,m 



With this choice we again arrive at formula (19), which is now applied to each 
component of the solution vector. 



2. 7 Roe 's approximate Riemann solver for systems 



Roe's approximate Riemann solver does not follow from working with relax- 



ation system (24), but instead from 



+ 



01 I 
-{Jf 2J 





f{q) - /i 



(27) 



In the above expression. 



dl 
dq 



where Q is the Roe average [23] and satisfies 

(g) iQr-Qi) = fiQr)-fm- 



dl 
dq 



(2J 



(29) 



With this choice the Riemann solution is obtained by splitting the jump be- 
tween the left and right states, {Qe, f{Qi)) and {Qr, f{Qr)), along the m dis- 
tinct eigenvectors of J in the following manner: 



Qr — Ql 
f{Qr) - fiQi 



a 



(30) 



where and are the p*^ eigenvalue and right eigenvector of the Roe matrix 
J, respectively. Just as in the scalar case, it seems as though the parameters a 
are overdetermined. However, since J satisfies the constraint (29), it can easily 



be shown that the first set of m equations involving Qr — Qi are identical to the 
second set of m equations involving f{Qr) — f{Qe)- In other words, there are 
only m distinct equations for m values of a; and therefore, a unique solution 
exists: 

aP = F-iQr-Qe), for p=l,...,m, (31) 
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where is the p^^ left eigenvector of J. This results in the following fluctua- 
tions: 



m 

1± 



^±AQ = ^ [s^f {P ■ (Q, - Q,)} r^. (32) 



Note that conservation follows from ( 29 ) . 



3 Residual distribution schemes 



We will describe in this section the basic residual distribution method for 
solving hyperbolic conservation laws in multidimensions. For further details 
we refer the reader to articles by Abgrall [1.2J and Abgrall and Mezine [lf5] . 
We consider a conservation law of the form 

+ V ■ /(g) = g,i + /(g),. + g{q),y = 0, (33) 

where {x, y) G are the spatial coordinates, t G M is the time coordinate, 
g G is the vector of conserved variable, and /(g) : — >■ M™^^ is the flux 
function. We will assume that this equation is hyperbolic, meaning that the 
m X m flux Jacobian matrix, 

'^H = ^-^^' (34) 

is diagonalizable with real eigenvalues for all n G such that = 1 and 
for all g in the solution domain. 



The first step in approximately solving (33) in some domain f2 C is to mesh 



the domain with a finite number of triangles. We will refer to this triangulation 
as Th, where h refers to a representative triangle radius, which in this work we 
just take to be the square root of the triangle area: h = \J\T\. Associated with 
this triangulation is a dual grid, which is constructed by connecting triangle 
centers to edge centers. A an example triangulation along with its dual grid 
is shown in Figure [1} 

Unlike the discontinuous Galerkin approach p^]|22] , approximate solutions are 
centered on triangle nodes (i.e., centers of the dual grid) rather than trian- 
gle centers. In order to obtain an update for these node centered values, we 
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integrate (33) over the median dual cell Ci and from t = t"^ to t = t""*"^: 



q{x,e+^)dx= 11 q{x,e)dx- f ff V-f{q)dxdt 
Ci JJc^ Jt" JJCi 

qix,f)dx- <f f(q)-dsdt. 

nTTl-n- Jt^ JdlCiDT) 



Next we define the median dual cell average and the time-averaged fluctuations 
through diCiHT): 

Q^^-L // q(x,e)dx, (35) 



Cj I J Jc.. 

$r = — / f /(g) ■ dsdt. (36) 
At Jt" JdidnT) 

Using these definitions, the update formula for a generic residual distribution 
scheme can be written as follows: 

Qr' = Q7-^,T. (37) 

In the remainder of this work, we will approximate the exact solution, q{x, t), 
with a piecewise constant representation, Q", that is constant on each median 
dual cell. We note that this view of RD schemes is slightly different than the 
standard view (e.g., P), where the approximate solution is usually viewed 
to be piecewise linear on each triangle T. Although these descriptions seem 
contradictory, in the case for first-order accuracy in time, both interpretations 
yield the same numerical schemes. The advantage of viewing the solution as 
being piecewise constant on each medial dual cell is that this naturally sets 
up a series of multidimensional Riemann problems in each triangle (see Figure 
[2]), which can be solved to construct the fluctuation^ In this way, we can 
then view approximate constructions of as approximate Riemann solvers. 

Computing the fluctuations $^ is generally done using the following frame- 
work (again, the two interpretations, piecewise constant on each dual cell vs. 
piecewise linear on each primal cell, make use of the exactly the same frame- 
work): 

(1) On each triangle T construct a total residual: 

$^ = // V-f^dx=(f / ^ ■ ds, (38) 



IT J&T 

where denotes an interpolant that passes through the three nodal 

In this work, the terms distributed residual and fluctuation mean the same thing 
and are used interchangeably. 
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values 

[xi, f{Qi)) for i = 1,2,3. 

For example if we simply use linear interpolation, the total residual can 
be written aJ^ 

= l^:fiQ^)■n^. (39) 
^ i=l 

Here represents the inward pointing normal vectors to the three edges 
of the triangle T. The length of Ui is equal to the length of the edge to 
which it is perpendicular. If the three nodes of triangle T are given by 
{xi, Hi) for z = 1, 2, 3, then the three scaled normals can be written as 



ni = 


(l/2 


- 1/3, 


Xs 


- X2 


n2 = 


(^3 


- 2/1, 


Xi 


- X3 


= 




- 2/2, 


X2 


— Xi 



(2) Once this total residual has been calculated, it is then distributed to each 
of the nodes of the triangle: 

The detailed strategy for how this distribution is accomplished yields a 
specific numerical method. 



3.1 Design principles for scalar conservation laws 



We first focus on design principles for scalar equations; in a subsequent sub- 
section we explain how to extend this to the systems case. In order to obtain a 



numerical update that produces a stable and accurate approximation to (33), 
we will need the distribution strategy to satisfy certain properties: 

(1) Numerical conservation: Since the interpolation of the numerical so- 
lution is continuous across triangle edges, conservation simply reduces to 
the following constraint: 

^<|.f = $^. (40) 

i=l 

In other words, in a given triangle, the sum of the distributed residuals 
must equal the total residual. 

(2) Monotonicity preserving: This condition makes sure that the numer- 
ical update satisfies a local maximum principle, which is needed to guar- 
antee that the update does not generate any new spurious maxima or 

^ For the standard N-scheme, which we will describe shortly, this is not the defini- 
tion of the total residual that is used. 
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minima. If we write 



<j^r = E4 Qr-Q- , (41) 



then the monotonicity requirement can be written as (see [13]): 

cjj>0 Vz,jG (1,2,3). (42) 



(3) Linear preserving: The order of accuracy of update (37) in the steady- 



state depends, among other things, on how accurately the total resid- 



ual (38) is calculated on each triangle [Ej. If we use formula (39), then 
= 0{h^) in the steady state, which is the correct order of accuracy if 
want an approximate solution, Qj, that is 0{h'^) accurate in the steady 
state. What we actually need in order to get an 0{h?) accurate steady 
state solution is that not only that = Oiji?), but that each distributed 
residual also satisfies $f = 0{h?). The distributed residuals can be writ- 
ten as 

$f = A$^ $^ = /?2$^ ^'I = /33^'^ (43) 
where /3j measures the fraction of the total residual that is distributed to 
node i. To ensure that the distributed residuals are of the same order as 
the total residual, we need to make sure that the /3j's remain bounded 
as /i ^ 0. Therefore, the 0{h'^) accuracy condition, or more commonly 
referred to as the linear preserving condition, can be writtten as follows 
(see 

/3i for ? = 1, 2, 3, is uniformly bounded independent of the mesh. (44) 



3.2 Scalar N-scheme 



Modern finite volume methods for hyperbolic PDEs are typically based on 
solving, either exactly or approximately, a Riemann problem between neigh- 
boring states. For multidimensional problems, a standard approach is to solve 
local ID Riemann problems and then use the information from the Riemann 
solutions to construct numerical fluxes or fluctuations (see Chapters 19-21 of 
LeVeque [20]). 

In the RD framework, however, a multidimensional Riemann problem is solved. 
In an arbitrary triangle T, we consider the Riemann problem between three 
constant states: Qi, Q2, and Q3 (see Figure |2]for a depiction). Exact solutions 
to multidimensional Riemann problems are at best expensive to evaluate, and 
in general not well-understood for many hyperbolic systems such as the Eu- 
ler equations from gas dynamics [2S]. Therefore, in practice an approximate 
method such as the N-scheme (the "N" stands for Narrow) [T^P^ is utilized; 
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this approach can be viewed as a multidimensional generalization of Roe's 
approximate Riemann solver 



Just as in the ID case, we define a Roe-average (henceforth called the Roe- 
Struijs-Deconinck average |12j): 



u 



dq 



(Q), 



(45) 



where Q is an average of the three nodal values Qi for i = 1,2,3 on the 
current triangle T. The Roe-Struijs-Deconinck average satisfies the following 



constraint, which generalizes the ID constraint given by (29): 

3 



<|.^ ^ V (^^ ■ n^ Q,= i f{q^) . ds, 



(46) 



where g'* is the linear interpolant passing through (xj, Qi) for i = 1,2, 3. If the 
fiux, /(g), is at most a quadratic function of q, then 



u 



r 



f'{Qi) + f'{Q2) + f'{Qz)) . 



(47) 



The approximate Riemann solution gives rise to the following set of fiuctua- 
tions: 



N-scheme: 



U ■ Uj 



{Qi - Q. 



(4J 



where is the so-called upwind parameter. In ID the the upwind param- 
eter relative to state Qi is always either Qi-i if m > or Qi+i if m < 0. In 



multidimensions, Qi, is obtained by enforcing the conservation constraint (|40|): 



3 



(49) 



where we have made use of the following two identities: 



3 3 
1=1 



= E [^'^ <5i + E [' 



i=l 
3 



■ fii 



i=l 



Qii 



u'^ ■ fii 



i=l 



+ E ■ 

1=1 



•T 



In order to demonstrate that the N-scheme is monotonicity preserving, we 



rewrite (48) in the form (41) with 



■ Hi 



> 0. 



(50) 
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From this it can be shown that update (37) is monotone under the following 
CFL condition: 



At < min 



2|a 



(51) 



3.3 Linear preserving limiters 



The N-scheme described so far is both conservative and monotonicity preserv- 
ing; however, it is not yet linear preserving. The problem with the previously 



described N-scheme is that the weights (43) are not uniformly bounded inde- 
pendent of the mesh. In order to modify the N-scheme to achieve uniformly 
bounded A's, Abgrall and Mezine [5] introduced a nonlinear limiting proce- 
dure. The limiting process takes the original A and replaces them with limited 
versions, denoted A. The simpler of the two approaches discussed in j|5| yields 
the following formulas: 

which guarantees that < A < 1. The limited residuals are then given by 

Limited N-scheme: <l>f = A (53) 

It is clear that this scheme is both linear preserving and conservative. Further- 
more, Abgrall and Mezine [5j proved that the limited N-scheme retains the 
monotonicity properties of the original N-scheme with the same CFL condition 
(|5T|. 



3.4 Extension to systems 



Following [12], the N-scheme is extended to systems of conservation laws by 
first defining the following averaged flux Jacobians: 

where Z is a parameterization of Q and 

Z = ^ (Zi + Z2 + Zs) . (55) 
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In order to achieve a conservative linearization we must find a parameter 
vector, Z, such that the following constraint is satisfied: 



E ■■T)Qr='f> Rz'') ■ ds, (56) 

where is the linear interpolant that passes through the points (xj, Zi) for 
i = 1, 2, 3 and 

- I (Z) Z,. (57) 



As was shown in [12] (see also p^fTT] ). constraint (56) will in general only be 
satisfied if we are able to find a parameterization, z = {z^{q), z^(g), . . . , z"^{q)), 
such that the flux, f{q{z)), depends at most quadratically on each z^ for 
p = 1,2, ... ,m. For the Euler equations from gas dynamics, as well as related 
systems, such a parameterization is known [T2] . 



Assuming that a parameterization has been found, we proceed by diagonaliz- 
ing the flux Jacobian: 

fii ■ J = RiKiRj^ ^, 

where Ri is the matrix of right eigenvectors and Aj is the diagonal matrix of 
eigenvalues. Following the philosophy of Roe's approximate Riemann solver 
[25] , the systems N-scheme is obtained by applying the scalar N-scheme to 
each characteristic component. This results in the following method: 

Systems N-scheme: = ^RiAfR;^ {Qi - Q^) . (58) 

The upwind parameter can be recovered by enforcing local conservation: 

-1 

' - " ' (59) 



Note that solving for involves inverting an m x m matrix. Finally, we 
note that although it is based on a generalization of the monotone scalar N- 
scheme, the systems N-scheme is in general only approximately non-oscillatory 
for nonlinear systems of conservation laws. In practice, however, this scheme 
has been shown to work quite well for steady-state shock computations for 
systems such as the Euler equations from gas dynamics [1]. 

The systems N-scheme described so far is not linear preserving. In order that 
the limiting procedure developed for scalar equations can be re-used for sys- 
tems, Abgrall and Mezine [5] proposed to project the distributed residuals into 
the eigenspace of the Roe-Struijs-Deconinck- averaged flux Jacobian in some 
direction n. In practice, the direction n is chosen from physical considerations. 
For example, in the case of the shallow water equations or the Euler equation 
from gas dynamics, an approach that gives good results in practice is to take 
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n to be the local Roe-Struijs-Deconinck-averaged fluid velocity: n = u. Once a 
limiting direction has been chosen, the limiting procedure can be summarized 
as follows: 



for p = 1 . . .m 

for i = 1,2,3: set = F ■ $f ; 



for i = 1,2,3 : set /3f 
for i = 1,2,3 : set /3f 



+ ' 



end 



for z = 1, 2, 3 : set = ^ /3f Gf r^, 

p=i 



where and are the p left and right eigenvectors of n ■ J, respectively. 



3.5 A correction for improved convergence 



As was pointed out by Abgrall [2], the N-scheme in conjunction with the 



limiting procedure outlined in Section 3.4 has one major drawback: the method 



does not in general converge to a steady-state solution. The problem is not 
with the N-scheme itself, since this method does converge to a steady-state, 
but instead the problem lies in how the N-scheme interacts with the limiting 
procedure. In the same paper, Abgrall [2] also provided a cure for this problem. 
He arrived at the following distributed residual: 

+e\T\-'^l^ Ki^'^ , (60) 

limited N-scheme correction 

where Ki = {ui ■ Tj /2, |T| is the area of triangle T, and ^ is a grid and 
solution dependent parameter. Notice that conservation is not affected by this 
correction term. In order to produce a scheme that converges to a steady-state 
solution, 6 needs to be chosen so that the correction is relatively small near 
shocks {6 = O (^|T|i/^^) and relatively large in smooth regions {0 = 1). Abgrall 
[2] proposed the following formula: 



^ = min 1, , ' , (61) 











+ 10-10 ; 



where ip^ is the projection of onto some important eigen-direction. In 
the case of the compressible Euler equations, ip^ should be taken to be the 
projection of onto the entropy wave. In Section [sj in which we consider 
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several numerical examples, we will refer to the limited and corrected N-scheme 
as the N\LC-scheme. 



4 Multidimensional relcLxation systems 



Having reviewed the relaxation scheme paradigm in Section [2] and residual 
distribution schemes in Section |3} we now turn to develop a multidimensional 
relaxation system framework. We again begin with the case of a scalar con- 
servation law and introduce the following relaxation system: 



(62) 





q 




q 

























1 

e 












,t 




,x 












where 

































1 







A^ = 










+ d 


1 















- d^; 








+ d') 




















1 






d^- 




d') 




\ + 




+ d') 






-c^d^ 











c^ + d^ 



(64) 



In these expressions, we assume that c,d& M^. Just as in the ID case, we 
will separate the effects of the hyperbolic left-hand side of this equation from 



(63) 



the relaxation source term on the right-hand side by by viewing (62) as being 
comprised of the following two sub-problems: 



(65) 











q 




+ A' 




+ A' 










,x 







1 




/_ 


e 


_9{q) - /^^_ 



(66) 



In subsequent discussion we will make use of the following matrix: 

A{n) = n^A^ + n^A^, 



(67) 
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where n G such that ||n|| 7^ 0. The three eigenvalues of A(n) are given by 

= ^ (c-n + cT- n) ± {n^d^ - n^c^f + {n^d^ - n^c^f , (68) 
A2 = 1 (c-n + cT-n) . (69) 



2 V / 2 

1 

2 

We note that these eigenvalues are strictly real. 



Next we define a relaxation procedure that is analogous to the ID case; the 
main difference is that in 2D a genuinely multidimensional Riemann problem 
such as the one depicted in Figure |2] must be solved. Instead of attempting to 
solve this exactly, we solve it with the standard N-scheme. The full procedure 
can then be summarized as follows: 

(1) Choose values for the parameters c^, c^, d^, and d"^. 

(2) On an arbitrary triangle, T, approximately solve the multidimensional 



Riemann problem associated with (65) by applying the standard N- 
scheme. 



(3) Approximate the effect of equation (66) on the solution calculated in 
Step (2) by directly setting fj,^ = f{q) and /i^ = g{q)- In other words, 
instantaneously relax the solution from Step (2) to the e limit. 



4-1 The N-scheme 



The first scheme that we will produce with the relaxation procedure is the 
N-scheme applied to the original scalar conservation law. We set 



c = d = u, 



(70) 



where u is the Roe-Struijs-Deconinck- average that satisfies (46). The above 
choice for c and d results in the following coefficient matrix for the relaxation 
system: 



A{n) 





-u^ {u ■ n) n^u^ + u ■ n 



n 



2 1 

n u 



-u [u ■ n) 
this matrix has eigenvalues given by 

A^'2'3 



1 2 
n u 



u ■ n, 



n^v? + u - n 



(71) 



(72) 



and, as in the ID case, has an incomplete set of eigenvectors (i.e., the eigen- 
values have algebraic multiplicity 3, but geometric multiplicity of only 2). 
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Because A{n) only has two linearly independent eigenvectors it cannot be diag- 
onalized; and instead, we reduce it to Jordan canonical form via the following 
similarity transformation: 



A{n) = SMS-\ 



(73) 



where 





n ■ u 


1 







1 


]_ 

n-u 





M = 





n ■ u 





and S = 




— 


— 










n ■ u 











(74) 



In approximately solving the Riemann problem via the N-scheme (Step (2) 
of the relaxation procedure), we will need to make sense of the expression 
[^(ri)]"'". Without a full set of eigenvectors, we do this in the following way: 



[Ain)]- 



SM-^S 



+ c-l 



where M 



[n ■ u\ 





1 



\n ■ u\ 



■ (75) 



This results in two possibilities: 



(1) [n-u] = {n ■ u) =^ 

(2) [n-M]+ = =^ [A{n)] 



[A{n)Y = A{n) 
n ■ u 

1 



n 



n 



u [n ■ u) 
{n ■ u) 



11 12 

u n u n 



2 1 2 2 

u n u n 



The first of these two expressions is exactly the result one should expect; the 
second expression, however, is somewhat troubling. We should expect that 
[Airii)]^ {Ui - f/^) = if {rii ■ m) < 0, where 

Ui = (Qi, f{Qi)) and = {Q^, pi^) . 

Instead, we are currently stuck with the following result when {n-u) < 0: 



[A(r?,)]+ (f/, -U,) = [-{u-n) {Q, -Q.)+n- {f{Qi) - /ll) } 



u 



u 



In order to clean up this result, we are forced to slightly modify the relaxation 



procedure for the N-scheme. We will leave the sub-problem (65) alone, but 



20 



replace sub-problem (66) with the the following system of ODEs: 





1 


/_ 


e 



(76) 



where g is a piecewise constant function in space that is constant on each 
triangle T; the value of this constant is Q, the multidimensional Roe-Struijs- 



Deconinck average (46). This results in the following modification of Step (3) 
in the relaxation procedure: 



(3) On each triangle T approximate the effect of (76) on the solution cal- 
culated in Step (2) by directly setting fi^ = u^q and /x^ = v?q. In other 
words, instantaneously relax the solution from Step (2) to the e ^ 
limit. 



Note that in general 
df 



dg 



q^{q)<1 7^ fiq) and (q) q ^ 9{q); 



and therefore, replacing (66) with (76) will yield a different numerical scheme. 



As we will demonstrate below, it is the scheme based on ( 76 ) that will repro- 
duce the N-scheme. 

The solution value at each node is now given by 

= (g„ u^Qi, u^Qi) . (77) 
Additionally, we enforce the condition: 

f/, = (Q„ n^Q„ u^Q,) . (78) 

With these modifications it is now true that [y4(?2j)]^ (Ui — f/^) = if (rij ■ u) < 
0. 

In order to proceed with the relaxation procedure, we solve a Riemann problem 



between three states of the form (77) with i = 1,2,3. Solving this Riemann 
problem via the N-scheme tells us that the residuals distributed to each node 
are given by 







otherwise. 



(79) 



where [/^ has the form (78). In order to determine [/^ in terms of the f/j values. 



we add the three residuals given by ( 79 ) and enforce that this sum yields the 
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total residual in triangle T: 



1=1 



i=l 



{Hi ■ u) Qi 
{Hi ■ u) Qi 
{Hi ■ u) Qi 



^0) 



^1) 



If > 0, then we note the following result on each triangle T: 

(1) 3 A; G (1, 2, 3) such that u ■ Hk > 0, 

(2) 3k e (1, 2, 3) such that u ■ Hk < 0. 

This result implies that there are two possibilities whenever ||n|| > 0: the 1- 
target case - 3 exactly one k s.t. u ■ Hk> 0, and the 2-target case - 3 exactly 
two k s.t. u ■ > 0. Without loss of generality, let us assume that u ■ Hi > 
and M ■ ns < 0, which yields one of the two possibilities: 

1- target solution: -u ■ ni > 0, m ■ n2 < 0, m ■ ns < 0, (82) 

2- target solution: u ■ ni > 0, m ■ ^2 > 0, u ■ < 0. (83) 



We consider each of these two cases below. 



4-. 1.1 The 1-target solution 

The 1-target case is easy to analyze: the total residual is completely distributed 
to the lone node that is downwind of the flow, which we have taken without loss 
of generality to be node 1. If we let $ denote the component of the residual 
corresponding to Q, then the 1-target case results in the following residual 
distribution: ^ 

which is the same result that one would obtain with the N-scheme on the 
original scalar conservation law. 



Jf..l.2 The 2-target solution 

The 2-target case involves distribution to two nodes, which we have taken 
without loss of generality to be the 1 and 2 nodes. From equation (81) we 
arrive at the following linear system that must be solved in order to obtain 
the upwind parameter f/^: 

- {A{Hi) + A{H2)) = ^(ns) U,. (85) 
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However, it is not difficult to show that 



A{ni) + A{n2) + ^(^3) = 



Ain^) = - (Ain,) + A{n2)) . (86) 



This imphes that f/^ = f/3. Therefore the 2-target case results in the following 
residual distribution corresponding to Q: 



^1 = 1 {ni ■ u) (Qi - Qs) , ^1 = ]^ {n, ■ u) (Qi - Q3) , $1 = 0. (87) 



This result is again identical to the original scalar conservation law. 



4-2 The RXN -scheme: genuinely multidimensional local Lax-Friedrichs 



One of the main difficulties with the N-scheme is that computing the upwind 
parameter for complicated systems of conservation laws can become pro- 
hibitively expensive. Despite this fact, few alternatives have been developed 
in the RD literature. One such alternative was introduced by Abgrall [2|, who 
considered a local Lax-Friedrichs-type method that was obtained, in analogy 
to the ID case, by taking the unstable "centered" residual and adding the ap- 
propriate numerical viscosity. In this work, we construct a new method based 
on the idea of relaxation systems; this scheme can be viewed as a different 
multidimensional generalization of the ID LLF method. For brevity we will 
call this method the RXN-scheme, which stands for "relaxation N-schem^^" 
In analogy with the ID LLF method as derived in Section 2.3, we make the 
following choice for the parameters c and d in (62)-(64): 



d 



Note that each triangle can have a different value of s ; this is why we call it 
a 'local' Lax-Friedrichs. This choice results in a coefficient matrix of the form 



A{n) 



n 



n 



nMs^) 



n2 fs^) 



^9) 



^ the words "N-scheme" appear here because we make use of the N-scheme to solve 
the homogeneous part of the relaxation system. 
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which has an eigenvector decomposition given by 



A{n) = SMS-^ 



-1 



Nil 



\\n\\ 



M 



\\n\\ 



\n\\s 



T 



\n\\s 



r 



2\\n\\s'^ 
„2 



2||n||s'^ 



2||n||s'^ 2||n|| 



(90) 



The Chapman-Enskog expansion for this relaxation system can be written to 
first order as 



q,t + fiq),x + giq),y ~ eV- 



{s^y-ifiq)? -f\q)9\q) 

-nq)9'{q) (s-^Y -{9'{q)f 



Vq . (91) 



The eigenvalues of the diffusion matrix in the above expression are 

A^ = (0', X' = {s^y-{f'{qr + 9'{qr), (92) 
which results in the following restriction on the choice of the lone parameter 

s^>\\f'{q)l (93) 



for all q E T. 



Applying the N-scheme to the relaxation system with coefficient matrix (89), 
yields the following residual 



where U^, = (Q^, ^l, /i^). Simplifying this expression gives 



(94) 



^js^Pill {Qi - Q*) + ■ {fiQi) - f^*) I 



(95) 



In order to calculate U^, we must enforce conservation: 



3 ]^ 3 



i=l 



i=l 



nt ■ f{Qt) 
n] is'')' 



{s'Y Q^ 



(96) 
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which results in the following linear system for the upwind parameters (Qj,, /i^): 



E 

i=l 



Til ^ 

s m,' 








1 2 



1 9 

n n 

2 1 







3 






= E 















n 



T 



n 



(4 ■ /(Q,) - s^Q^ 



\nj\\Qj 

n 
\\n 



(97) 



The solution to this linear system can be written as 



-yy 

1 3 3 

EE 



N ^ 



where 



i=i i=i 



N 



\ni\ 



■fiQj)-s'Q. 



n]n\ 



EE 

i=i j=i 



1122 1122' 

nl n- Hj rij — nl 



\ni\ 



(98) 
(99) 

(100) 
(101) 



Let us now take a moment to reflect on what just happened. Although the 
original coefficient matrix, (89), for this method was comically simple, after 
applying the N-scheme to this system on an arbitrary triangle, the resulting 
upwind parameters are somewhat complicated. On the other hand, we see 
from equation (97) that the parameter is completely decoupled from /Ij,. 
We make use of this last fact to construct an alternative scheme in the following 
way: instead of computing the components of /I^ from (99)-(101), we enforce 

= RQ.) (102) 

by again invoking the e — > relaxation limit. Although this direct enforcement 
clearly gives a different scheme than if we had used (99 )-( 101 ), what we achieve 
with this approach is a very simple method that we refer to as the RXN-scheme 
(relaxation N-scheme). In terms of the residual distributed to node i in the 
Q-variable, we now obtain the following expression: 

RXN-scheme: = \s^ \\ni\\ {Qi -Q.) + \ni- [f{Qi) - f{Q,)) , (103) 



where is given by (98). Note that this method is automatically conservative 
since still satisfies the first equation in linear system (97). 

Theorem 4.1 (Monotonicity) If there exists a Q such that 

3 3 

E ■ fiQ^) = E ■ fiQ) (104) 



i=l 
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and 



> max 



{iimii, ii/'(Q2)ii, iim)ii}, (105) 

where for each j 

n, ■ [Kq,) - KqS) = ■ r (q,) (g, - g.) (loe) 

from the Rankine-Hugoniot conditions, then the RXN-scheme as defined by 
(103) and (98) satisfies the following condition: 



r 



E {Q^ - Q 



JJ ' 



(107) 



where > for all i,j = 1, 2, 3. This condition along with the following CFL 
constraint on the time step: 

At < mill 



2 1 a 



(108) 



Ell ^ II T f ' 
T-.ieT iFill s ) 

is enough to guarantee that the RXN-scheme with forward Euler time-stepping 
(31) is monotonicity preserving. 



Proof. (1) Using the Rankine-Hugoniot conditions ( 106 ), we rewrite the RXN- 
scheme as 



= (Qi-Q,), where P, 



Til ^ 
S Hi 



n, ■ / ' (g. 



(109) 



Similarly, we rewrite (|98|) as follows 

g 



V^"ijr' ^ where NJ ^ s^\\n,\\ - ■ f\Q) 



;iio) 



Note that the above expression was obtained by makin g use of ( 104 ) and the 
identity: Y%=i^k ■ f'iQ) = 0. Combining expressions (109) and (110) yields 
(irOTl) with 



We note that ^ > Vi,j G (1,2,3), because (105) implies that Pj > Q 
Vi G (1, 2, 3) and Nj > V j G (1, 2, 3). 



(2) We now insert expression (107) into (37) and simplify: 

At 



gr' = g 
gr' = 



1 



\C^\ 

At 



E E4(Qr-g-), 
E E4|Qr + ^ E E4Q"- 



T-.ieT jer 
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Monotonicity is achieved if Q^^^ is a convex average of all of the surrounding 



Q^. Since each > 0, we obtain a convex average provided that 




< 1 



At, < 



\Ci 



At < min 



J- 



Er Ej jv^T 



E 



The time restriction is clearly satisfied if we take (108). | 

In practice the time step presented in the above theorem is overly restrictive. 
In the numerical simulations presented in Section |5j we instead use the same 
time-step as used with the N-scheme: 85% of the maximum CFL number given 
by expression (51). 



4-3 Systems N-scheme 



The systems generalization of coefficient matrix (89) for a system of m con- 
served variables is the following 3m x 3m matrix: 



A{n) 



01 



n 



n 



.J).P 



n J + n ■ J 



n^.P + n ■ J 



;ii2) 



where I is again the mxm identity matrix, 01 is the mx m matrix with zeros 
in every entry, and J = (^J^, J^j is the flux Jacobian matrix evaluated at the 



Roe-Struijs-Deconinck average [12]. The systems generalization of (77)-(78) 
are the following 3m x 1 vectors: 



Ui = (Qi, pQi, PQ,) and = (Q., J^Q., PQ, 



113) 



In order to calculate the appropriate residuals in the relaxation procedure, we 
need to again understand how to create the matrices [^(n)]''' and [A(n)]~. As 
in the scalar case this is complicated by the fact that A{n) does not have a 
full set of eigenvectors. In particular, the Jordan canonical form of this matrix 
can be written as 



A{n) = S 







xp 


1 







S-\ where = 





XP 















XP 



;ii4) 
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Here is the p^^ eigenvalue of the mxm matrix n- J. We omit the comphcated 
expression for the matrix S. In order to obtain an expression for [A(n)]''', one 
has to replace each in the above expression with [A^]^. Carrying this out 
results in the following matrix: 



[A{n)r 



-J- 



J\j- 
.PJ- 



v}P + P 
n^P 



n^P 
n^P + P 



(115) 



where = (ji ■ Pj = RA^R ^, A is the diagonal matrix of eigenvalues 

of n ■ J, and R is the corresponding matrix right-eigenvectors. An analogous 
formula for [y4(n)]~ can also be readily constructed. From the above expression 
we find that 



[A{n)f U, 



(n-J^Q,, P{n-j)^Qi, P (n ■ j)^ Q, 



1 1 



(116) 



Having established expressions for [y4(fi)]^, we now proceed by applying the 
N-scheme to the relaxation system: 



i=l 



where 



j:[Mnr] u., (117) 



E [Mm 



i=l 



Si Ji 
Ji J ~^ J Ji 



\ i=l 



Y.^J- 

Y.^J^ 



:ii8) 



and Ji = (jii ■ Pj. The unique solution to the linear system in (117) is Ui, 

(q., Pq,, Pq,) with 



Q^ = \T.{^i-J) ) {Y.{^^-J) Q^)^ (119) 

and the component of the residual yjj associated with Q can be written as 

= \{n,-Pi^ {Q,-Q,). (120) 
This result shows that this relaxation scheme identically reproduces the sys- 



tems N-scheme (58)-(59). 
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4-4 Systems RXN-scheme 



Just as the LLF method in the one-dimensional case, the RXN-scheme extends 
to systems of conservation laws in a simple manner. All that we have to do 
is apply the scalar version of the scheme to each component of the vector 
conserved variable. The coefficient matrix in the relaxation procedure can be 
written as 



'1211 



where I is the mxm identity matrix. In order to satisfy the sub-characteristic 
condition we require that 





01 n^I 




A{n) = 


n^s'^I 01 


01 




n^s^I 01 


01 



s> max \J[XP'''{q)] +[XP'y{q)y 

p=l,...,m 



(122) 



over all g e T. In the above expression A^'^ and A^'^ are the p^^ eigenvalue of 
df/dq and dg/dq, respectively. 

We note that the systems RXN-scheme, and in particular, the version of this 



scheme with limiters (Section 3.4) and convergence corrections (Section 3.5), 
provides an alternative to the systems N-scheme that does not require the 
inversion of an m x m matrix in each element at each time level, nor does it 
require any special entropy fixes or special treatment near stagnation points. 
Since this method is also simpler than the N-scheme, it should also yield some 
gains in computational efficiency. The systems N-scheme and RXN-scheme are 
compared in detail in Section [5j 



4-5 RXN-scheme in d-dimensions 



The above procedure for obtaining the 2D RXN-scheme can be generalized 
to any space dimension. In the ci-dimensional case we arrive at the following 
scheme: 

RXNd -scheme: <I>f = — ^ ^ — ^, (123) 

where 

In particular, we note that for c? = 1, this scheme exactly reduces to the ID 
local Lax-Friedrichs method ^24j . We also mote that the 1/d geometric factor 
comes from the d-dimensional N-scheme; see for example equation (7) in jTT] . 



= ' ' w;;.;„^n • (124) 
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5 Numerical examples 



In this section we compare the N-scheme and the newly proposed RXN scheme 
on several numerical examples. We will refer to the versions of the N-scheme 
and RXN-scheme that have been limited according to the procedure outlined 
in Section 3l4] and corrected according to the procedure outlined in Section 



3.5 as the N\LC-scheme and RXN\LC-scheme, respectively. 



5.1 Steady-state advection 

First, we consider the advection equation on [0, 1] x [0, 1]: 

g,f + M-Vg = 0, (125) 
with non-divergent velocity and boundary conditions given by 

u{x,y) = {-ny, tcx) , 

g(l,y)=0, g(x,0) 



sin^TT (^)) if 0.1 < X < 0.7, 
otherwise. 



This same problem was considered in [2]. 

For a non-divergent velocity field, an elegant way to solve the advection equa- 
tion using the N-scheme is through the introduction of a streamfunction: 



2 

such that u = {dip/dy, —dip/dx). The N-scheme can then be written as 

*r = H {Qi - Q.) , 

where 

^1 = ^ (^(3^2, 1/2) - i'i.x^, 2/3)) , 

k2 = ^{ip{x3,y3) -tp{xi,yi)), 
h = ^{i^{xi,yi) -i){x2,y2))- 

The advantage of this formulation is that we achieve numerical conservation, 
even though the equations are solved in advective form. 
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For the RXN scheme, we use residual (103) where the flux functions are given 
by 



f{Qi) = UiQi, t[Q*) = u^Q*^ ^ II : 

Pi + P2 + P3 



and Qi, is given by (97) as usual. 



Results on a grid with 5592 elements and 2903 nodes is shown in Figure |3] 
displayed in each panel are (a) the basic N-scheme, (b) the basic RXN-scheme, 
(c) the limited N-scheme (no convergence correction is needed for the limited 
N-scheme on scalar equations), and (d) the RXN\LC-scheme (convergence 
corrections are needed for the limited RXN-scheme, even for scalar problems). 
These results show that the basic RXN scheme is far more diffusive than the 
N-scheme. However, with limiting and convergence corrections, the RXN\LC- 
scheme gives results comparable to the limited N-scheme. Convergence histo- 
ries for the limited N-scheme, limited RXN-scheme, and the RXN\LC-scheme 
are shown in Figure |4j 

Clearly there is no advantage in using the RXN\LC-scheme over the limited 
N-scheme for a scalar problem, since the two methods have the same com- 
putational cost and the scalar limited N-scheme does not require convergence 
corrections. However, for hyperbolic systems such as the Euler equations, the 
RXN\LC-scheme provides a simpler algorithm with lower computational cost 
than the N\LC-scheme. 



5.2 Transonic flow around the NACA 0012 airfoil 



Next we consider transonic flow around the NACA 0012 airfoil using the com- 
pressible Euler equations from gas dynamics as our model. The Euler equations 
can be written as 









p 




pu 


pu 


+ V- 


puu + pi 


£ 




u{£ +p) 



where p is the fluid density, u = (u^, m^) is the fluid velocity, £ is the total 
energy, and p is the fluid pressure. The ideal gas law closes the system by 
relating the pressure to the other fluid variables: 

£ = ^ + lp\\uf, (127) 
7 — i / 

where 7 is the ideal gas constant. In this example we take 7 = 1.4. 
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The computational domain and numerical grid is shown in Figure [Hj^a). In 
Figure |5](b), we show a zoomed- in view of the numerical grid near the airfoil. 
The boundary conditions are such that subsonic flow with Mach number 0.85 
enters from the left boundary at an angle of +1° from the horizontal axis. As 
the flow impinges on the airfoil, two supersonic bubbles are created above and 
below the airfoil. The supersonic flow is decelerated to the ambient subsonic 
flow through the creation of two shock waves (again, one above the airfoil 
and one below) . This problem has been considered in several papers including 
IT1I2]. 



Isolines of the Mach number and the pressure are shown in Figure [6j We also 
plot in Figure [7] the Mach number along the top and bottom edges of the 
airfoil. From this figure we note that the location of the shocks is in very good 
agreement between the two methods, while the location at the front of the 
airfoil where the solution goes from subsonic to supersonic is slightly different 
for each method. Furthermore, the RXN\LC-scheme is more diffusive than 
the N\LC-scheme, which results in slightly more entropy production near the 
airfoil for RXN\LC than N\LC. This can be seen both in Figure [6| where 
we see bending of the Mach isolines near the airfoil, as well as in Figure |8j 
Overall, however, both of these figures indicate remarkable agreement between 
the two solutions. In particular, the RXN\LC solution is far closer to the N\LC 
solution than the MUSCL-type scheme that was presented in [1]. We also note 
that the RXN\LC-scheme runs twice as fast as the N\LC-scheme. 

Finally, in Figure [9] we show the L2-norm of the total residual as a function of 



time. We note that the fix of Abgrall [2] (see Section 3.5 ) is critically important 
in bringing both methods to convergence. Without this fix both methods stall 
at a total residual of only about 10~^. We also find that the RXN\LC-scheme 
has a slightly better convergence rate than the N\LC-scheme. 



5. 3 Supersonic flow around a cylinder 



Next we consider flow around a cylinder with Mach number M.00 = 5. The 
computational domain is [—2, 0] x [—3, 3] with a cylinder of unit radius centered 
at (0, 0). In this example we found that we needed to run all the schemes at a 
CFL number of 0.4 in order to obtain well-converged results. The steady-state 



pressure on a grid with 5144 elements and 2656 nodes is shown in Figure 10 



for the following schemes: (a) the basic N-scheme, (b) the basic RXN-scheme, 
(c) the N\LC-scheme, and (d) the RXN\LC-scheme. The basic RXN-scheme 
solution is far more diffusive than the basic N-scheme, but once limiters and the 
convergence corrections are added, the N\LC and RXN\LC schemes produce 
comparable results. In fact, the RXN\LC-scheme converges faster than the 



N\LC-scheme, as can be seen in the convergence plot in Figure 11 Finally, we 
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note that the RXN\LC-scheme again runs about twice as fast as the N\LC- 
scheme. 



5.4 Subsonic flow around a cylinder 



Finally, we consider flow around a cylinder with Mach number A4oo = 0.35. 
This problem has been considered in several papers including The com- 
putational domain is [—7,7] x [—7,7] with a cylinder of radius 1/2 centered 
at (0,0). The steady-state Mach number on a grid with 12552 elements and 



6404 nodes is shown in Figure 12 for the (a) N\LC and (b) RXN\LC schemes. 
Near the cylinder both methods produce comparable results. Away from the 
cylinder the grid resolution becomes coarser; and therefore, visible differences 
in the two methods appear. In these regions the RXN\LC scheme produces 
slightly more diffused contours than the N\LC scheme. 



Shown in Figure 13 are the deviation of the physical entropy, s = \og{p/p'^), 
from the ambient entropy, Soo = log(l/7'^): S = (s — Soo)/|soo|- Panel (a) is 
the N\LC-scheme and panel (b) is the RXN\LC-scheme. The minimum and 
maximum values of E for the N\LC and the RXN\LC schemes are (—4.101 x 
10"^ 4.324 X 10"^) and (-1.471 x 10"^ 1.291 x 10"^), respectively. Each panel 
consists of 31 contour lines ranging from the minimum to the maximum S for 
each scheme. Therefore, these results show that the RXN\LC-scheme has a 
smaller entropy deviations, but that this error is more spread out behind the 
cylinder, while the N\LC-scheme has larger entropy deviations, but that this 
error is more concentrated near the x-axis. 

The total residual as a function of time is shown in Figure [T4j Both methods 
give essentially the same convergence rates for this example. Finally, we note 
that the RXN\LC-scheme again runs about twice as fast as the N\LC-scheme. 



6 Conclusions 



In this work we have extended the results of LeVeque and Pelanti [21] to gen- 
uinely multidimensional residual distribution schemes. Speciflcally, we have 
shown that the N-scheme, both the scalar and the systems version, can be 
derived from a relaxation principle. Furthermore, using a genuinely multidi- 
mensional extension of the ID local Lax-Friedrichs relaxation principle, we 
have derived a novel residual distribution scheme. The main beneflt of this 
approach is that it does not require the inversion of an m x m matrix, where 
m is the number of conserved variables, at each time step in each grid ele- 
ment. The new method also does not require the use of Roe-Struijs-Deconinck 
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averages. Using several examples of the 2D Euler equations from gas dynam- 
ics, including an example of transonic flow around the NACA 0012 airfoil, 
supersonic flow around a cylinder, and subsonic flow around a cyhnder, we 
have compared the limited and corrected N-scheme (N\LC) with the newly 
proposed scheme (RXN\LC). These comparisons show that despite being com- 
putationally less expensive, the new method is capable of producing results 
comparable to those of the N\LC-scheme, although often with slightly more 
numerical diffusion. For more complicated equations such as magnetohydro- 
dynamics or the general relativistic Einstein equations, we believe that the 
benefit of a simpler and computationally less expensive algorithm will far out- 
weigh the slight increase in numerical dissipation. We will consider some of 
these more complicated systems in future work. 

Finally, we would like to point out that the numerical code used in this work, 
including all of the numerical grids, will be made publicly available as part of 
the REDPACK software project. For more information see 

http : //www . math . wise . edu/ ~rossmani/sof tware . html . 



Acknowledgements. The author would like to thank the two anonymous 
reviewers for their very helpful comments. This work was supported in part 
by NSF grants DMS-0619037 and DMS-0711885. 



References 

[1] R. Abgrall. Toward the ultimate conservative scheme: Following the quest. J. 
Comp. Phys., 167:277-315, 2001. 

[2] R. Abgrall. Essentially non-oscillatory residual distribution schemes for 
hyperbohc problems. J. Comp. Phys., 214:773-808, 2006. 

[3] R. Abgrall. Residual distribution schemes: Current status and future trends. 
Comput. Fluids, 35:641-669, 2006. 

[4] R. Abgrall and M. Mezine. Construction of second order accurate monotone and 
stable residual distribution schemes for unsteady problems. J. Comp. Phys., 
188:16-55, 2003. 

[5] R. Abgrall and M. Mezine. Construction of second-order accurate monotone 
and stable residual distribution schemes for steady problems. J. Comp. Phys., 
195:474-507, 2004. 

[6] R. Abgrall and P.L. Roe. High order fluctuation schemes on triangular meshes. 
J. Sci. Comp., 19:3-36, 2003. 



34 



[7] F. Bouchut. Nonlinear Stability of Finite Volume Methods for Hyperbolic 
Conservation Laws and Well-Balanced Schemes for Sources. Birkhauser Verlag, 
2005. 

[8] G.Q. Chen, CD. Levermore, and T.P. Liu. Hyperbolic conservation-laws with 
stiff relaxation terms and entropy. Comm. Pure Appl. Math., 47:787 830, 1994. 

[9] G. Chesshire and W.D. Henshaw. Composite overlapping meshes for the 
solution of partial-differential equations. J. Comp. Phys., 90:1-64, 1990. 

[10] B. Cockburn and C.-W. Shu. The Runge-Kutta discontinuous Galerkin method 
for conservation laws V: Multidimensional systems. J. Comp. Phys., 141:199- 
224, 1998. 

[11] A. Csfk, M. Ricchiuto, and H. Deconinck. A conservative formulation of the 
multidimensional upwind residual distribution schemes for general nonlinear 
conservation laws. J. Comp. Phys., 179:286-312, 2002. 

[12] H. Deconinck, P.L. Roe, and R. Struijs. A multidimensional generalization of 
Roe's difference splitter for the Euler equations. Comput. Fluids, 22:215-222, 
1993. 

[13] H. Deconinck, G. Struijs, G. Bourgeois, and P. Roe. Compact advcction schemes 
on unstructured meshes. In VKI Lecture Series 1993-04, Computational Fluid 
Dynamics, 1993. 

[14] A. Harten, P.D. Lax, and B. van Leer. On upstream differencing and godunov- 
type schemes for hyperbolic conservation laws. SIAM Review, 25:35-61, 1983. 

[15] C. Helzel. A high-resolution rotated grid method for conservation laws with 
embedded geometries. SIAM J. Sci. Comput, 26:785-809, 2005. 

[16] G.S. Jiang and C.-W. Shu. Efficient implementation of weighted eno schemes. 
J. Comp. Phys., 126:202-228, 1996. 

[17] S. Jin and Z.P. Xin. Relaxation schemes for systems of conservation-laws in 
arbitrary space dimensions. Comm. Pure Appl. Math., 48:235-276, 1995. 

[18] A. Kurganov and E. Tadmor. New high-resolution central schemes for nonlinear 
conservation laws and convection-diffusion equations. J. Comp. Phys., 160:241- 
282, 2000. 

[19] R.J. LcVcquc. Wave propagation algorithms for multi-dimensional hyperbolic 
systems. J. Comp. Phys., 131:327-335, 1997. 

[20] R.J. LcVcquc. Finite Volume Methods for Hyperbolic Problems. Cambridge 

University Press, 2002. 

[21] R.J. LcVeque and M. Pelanti. A class of approximate Riemann solvers and their 
relation to relaxation schemes. J. Comp. Phys., 172:572-591, 2001. 

[22] W.H. Reed and T.R. Hill. Triangular mesh methods for the neutron transport 
equation. Technical Report LA-UR- 73-479, Los Alamos Scientific Laboratory, 
1973. 



35 



[23] P.L. Roe. Approximate Riemann solvers, parameter vectors, and difference 
schemes. J. Comp. Phys., 43:357-372, 1981. 

[24] V.V. Rusanov. Calculation of interaction of non-steady shock waves with 
obstacles. J. Comp. Math. Phys. USSR, 1:267-279, 1961. 

[25] C.W. Schulzrinne, J.P. Collins, and H.M. Glaz. Numerical-solution of the 
Riemann problem for 2-dimensional gas-dynamics. SIAM J. Sci. Comput., 

14:1394-1414, 1993. 

[26] E. van der Weide and H. Deconinck. Positive matrix distribution schemes for 
hyperbolic systems. In Computational Fluid Dynamics 1996, pages 747 — 753. 
Wiley, 1996. 

[27] B. van Leer. Towards the ultimate conservative difference scheme V. A second 
order sequel to Godunov's method. J. Comp. Phys., 32:101-136, 1979. 



36 



O Triangle center 

# Triangle node 

• Center of triangle edge 



Grid lines of median dual cell; 

Grid lines 



Fig. 1. Sample triangulation and dual grid. 




Fig. 2. A depiction of the multidimensional Riemann problem that must be solved 
in each triangle. The numerical solution is piecewise constant on each median dual 
cell. For example, the approximate solution on the three dual cells that overlap 
the triangle shown in this figure are Qi, Q2, and Q3. Note that the area of each 
of the three sections is the same, the midpoint where the dashed lines meet is 
(xi + X2 + X3) /3, fik are the inward-pointing normal vectors to each edge, and the 
magnitude of nfc is equal to the length of the edge to which it is orthogonal. 
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Fig. 3. Advection equation example. Shown in these panels are (a) the basic 
N-scheme, (b) the basic RXN-scheme, (c) the limited N-scheme (no convergence 
correction is needed for the limited N-scheme on scalar equations), and (d) the 
RXN\LC-scheme (convergence corrections are needed for the limited RXN-scheme, 
even for scalar problems). These results show that the basic RXN scheme is far more 
diffusive than the N-scheme. However, with limiting and convergence corrections, 
the RXN\LC gives results comparable to the limited N-scheme. 
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L2-norm of residual vs. time 




Fig. 4. Advection equation example. L2-norms of the total residual for the limited 
N-scheme, the limited RXN-scheme, and the RXN\LC-scheme. This plot shows that 
the limited N-scheme does not need additional convergence corrections for scalar 
equations, but the limited RXN-scheme clearly does. 



Numerical Grid Numerical Grid 




Fig. 5. The numerical grid for the NACA 0012 example. Panel (a) shows the entire 
domain, while Panel (b) shows a zoomed-in view of the airfoil. The grid has a total 
of 14,284 elements and 7,298 nodes. The smallest grid elements near the airfoil have 
a triangle radius h that is roughly 30 times smaller than that of the largest grid 
elements on the outer boundaries. 
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Mach number Mach number 




(a) -1 1 2 (b) -1 1 2 



Pressure Pressure 




(c) -1 1 2 (d) -1 1 2 

Fig. 6. Steady-state solution of transonic flow around the NACA 0012 airfoil. 
Panels (a) and (b) show isolines of the Mach number for the N\LC-scheme and 
RXN\LC-scheme, respectively. Panels (c) and (d) show isolines of the pressure for 
the N\LC-scheme and RXN\LC-scheme, respectively. 
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Fig. 7. The Mach number along the top and bottom edges of the NACA 0012 airfoil. 
Panel (a) is the N\LC-scheme solution, while Panel (b) is the RXN\LC-scheme 
solution. In Panels (c) and (d) we directly compare the Mach number profiles for 
each method: Panel (c) shows the Mach number on the top portion of the airfoil and 
Panel (d) shows the Mach number on the bottom portion of the airfoil. These results 
show that the RXN\LC-scheme produces comparable results to the N\LC-scheme. 
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1- in o 113 1- 

(a) do. 

Fig. 8. Deviation of the physical entropy, s = log{p/ p'^), from the ambient entropy, 
Soo = log(l/7'''): S = (s — Soo)/|soo|- Panel (a) is the N\LC-scheme and panel (b) is 
the RXN\LC-scheme. These results again show that the RXN\LC scheme is slightly 
less accurate than the N\LC scheme near the airfoil. The same contour values are 
plotted in each panel: S = 0.002 : 0.002 : 0.08. The minimum and maximum values 
of S for the N\LC and the RXN\LC schemes are (-4.906 x lO""^, 8.029 x lO^^) and 
(-2.778 X 10-^ 5.805 x lO'^), respectively. 
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L2-norm of residual vs. time 
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Fig. 9. L2-iiorm of the total residual as a function of time for the NACA 0012 



problem. We note that the fix of Abgrall [2] (see Section 3.5 ) is critically important 
in bringing both methods to convergence. Without this fix both methods stall at a 
total residual of only about 10~^. 
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Fig. 10. Supersonic flow past a cylinder with Moo = 5. Shown in these panels 
are the steady-state solutions as computed with the (a) basic N-scheme, (b) basic 
RXN-scheme, (c) N\LC-scheme, and (d) RXN\LC-scheme. These results show that 
the RXN solution is much more diffusive than the N-scheme solution; however, 
once the limiters and convergence corrections are included, the N\LC and RXN\LC 
schemes produce comparable results. 
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Mach number 




Fig. 12. Subsonic flow past a cylinder with A4oo = 0.35. Shown in these panels are 
isolines of the Mach number for the (a) N\LC and (b) RXN\LC schemes. Near 
the cylinder both methods produce comparable results. Away from the cylinder the 
grid resolution becomes coarser; and therefore, visible differences in the two meth- 
ods appear. In these regions the RXN\LC scheme produces slightly more diffused 
contours than the N\LC scheme. 
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Deviation of physical entropy 
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Fig. 13. Subsonic flow past a cylinder with Moo = 0.35. Shown in these pan- 
els are the deviation of the physical entropy, s = \og{p/ p^), from the ambient 
entropy, Sqo = log(l/7'^): S = (s — Soo)/|soo|- Panel (a) is the N\LC-scheme 
and panel (b) is the RXN\LC-scheme. The minimum and maximum values of S 
for the N\LC and the RXN\LC schemes are (-4.101 x 10"^ 4.324 x 10"^) and 
(—1.471 X 10~^, 1.291 X 10~^), respectively. Each panel consists of 31 contour lines 
ranging from the minimum to the maximum S for each scheme. Therefore, these 
results show that the RXN\LC-scheme has a smaller entropy deviations, but that 
this error is more spread out behind the cylinder, while the N\LC-scheme has larger 
entropy deviations, but that this error is more concentrated near the x-axis. 



L2-norm of residual vs. time 
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Fig. 14. L2-norm of the total residual as a function of time for subsonic flow past a 
cylinder. Both methods give essentially the same convergence rates for this example. 
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